1 Goal

Now that we have a probably working method, we need to start validating the generated ensemble of Tgav time series vs the original ensemble. As Claudia has pointed out, there’s two elements of this:

1. We have to make sure that the jumps between pasted windows aren’t too … well ‘wrong’ for some definition of wrong. We do also need to think about how the variability within each pasted window compares across pasted windows sort of, especially if they’re coming from different RCPs. For annual, global temperature, it is probably fine but that’s something we may have to think about more carefully when we look at the monthly gridded netcdfs (I think it would show up in the power spectra of eofs maybe?).

2. We have to do the standard comparison of training ensemble to generated ensemble along the lines of Tebaldi et al 2020.

The reality is that because this is a top down constructive method, validation is so much more of a vital step than it was with, say, fldgen. Fldgen is mathematically guaranteed to produce new realizations that will validate well and the validation is just a formality once the emulation and generation is complete. For stitches, it’s an absolutely integral part of the generation step.

We will start with 2, as it’s an implicit partial validation of 1. Then we will look in more detail at 1. This notebook will be concerned with more traditional methods of validation.

We will also be considering preliminary exploration of 2 hyper parameters: whether the target experiment is included in the archive at all, and the neighborhood size (tol) for matching. Because we want to think about the errors of our method (and how those relate to the hyperparameters) and not just a single draw of a collapse free generated ensemble like might be used in practice, we will be doing many many draws of max-sized collapse free generated ensembles so that we can aggregate the errors across those draws. In addition to looking at errors seen in a single draw.

2 Set Up

set.seed(1)
# required packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(knitr)

# Then load the functions we will want to use, these are currently written in R and will be translated into python. 
source("functions/nearest_neighbor_matching.R") # the function match_neighborhood is defined here
source("functions/permuting_matches.R")  # the function permute_stitching_recipe is defined here
source("functions/stitching_functions.R")  # the function stitch_global_mean is definend here 
source("functions/gridded_recipe.R") # wrappers for the work flow from matching to output csv needed by python for 
                                     # layering in gridded data to output as netcdf

# How many times do we want to draw full generated ensembles:
Ndraws <- 50

Load the inputs that will be used.

# Time series of the raw global mean temperature anomaly, this is the data that is 
# going to be stitched together. Historical data is pasted into every scenario. No
# smoothing has been done.
tgav_data <- read.csv("inputs/main_raw_pasted_tgav_anomaly_all_pangeo_list_models.csv", 
                      stringsAsFactors = FALSE) %>% dplyr::select(-X)


# A chunked smoothed tgav anomaly archive of data. This was tgav_data but several 
# steps were taken in python to transform it into the "chunked" data and save it 
# so that we do not have to repeate this process so many times. 
archive_data <- read.csv('inputs/archive_data.csv', stringsAsFactors = FALSE)


# The target data - all ssp245 realizations (smoothed because it's the target for
# matching)
archive_data %>%
  filter(experiment == 'ssp245') ->
  target_data


# The corresponding unsmoothed original data for validation
tgav_data %>%  
  filter(model == unique(target_data$model) & experiment == unique(target_data$experiment))  %>%
  select(-file, -timestep, -grid_type) -> 
  original_data

3 Generate the new ensembles

We’re actually going to generate multiple different ensembles for a target of ssp245, based on two different archives.

Archive 1 will exclude the ssp245 runs from potential matches.

Archive 2 will include the ssp245 runs - when we get to the point of targeting SSP585 runs, I think we will have to to include those runs in the archive to get at all decent generated Tgavs at the end of the century when temperature is higher. It’s still not just regurgitating the target ensembles, it will include mixing and matching from those parts.

We will also do some light exploration of the matching neighborhood size (tol argument).

In addition to SSP534-over having enough data issues to fully exclude, some of the realizations for SSP434 and SSP460 just don’t have future data at all. We need to go back to pangeo and fix this for sure, but for now, I’ll work with the archive I have. Because we don’t want to lose too much of the archive now, I will only remove those realizations that don’t have all the years of data.

archive_data %>% 
  # Only match to the same ESM:
  dplyr::filter(model == unique(target_data$model)) %>%  
  # Don't match to any ensemble members in the target experiment:
  dplyr::filter(experiment != unique(target_data$experiment) ) %>%
  # Exclude ssp534-over:
  dplyr::filter(experiment != 'ssp534-over') %>%
  # eliminate any individual realizationXexperiment runs that don't have all the years:
  group_by(experiment, variable,   model  , ensemble) %>%
  mutate(nWindows = n()) %>%
  ungroup %>%
  filter(nWindows == 28) %>%
  select(-nWindows) ->
  archive_subset1 


archive_data %>% 
  # Only match to the same ESM:
  dplyr::filter(model == unique(target_data$model)) %>%  
  # Exclude ssp534-over:
  dplyr::filter(experiment != 'ssp534-over') %>%
  # eliminate any individual realizationXexperiment runs that don't have all the years:
  group_by(experiment, variable,   model  , ensemble) %>%
  mutate(nWindows = n()) %>%
  ungroup %>%
  filter(nWindows == 28) %>%
  select(-nWindows) ->
  archive_subset2 

** TODO ** look into NA’s getting stitched in but for now, just drop any generated ensemble member that has any NA values; probably related to the order of matching in and discarding archive points since this analysis is pushing the limits of the functions we wrote a bit further than the stress tests we did while writing them. Some might were related to the 2014 as a window year in SSP434 and SSP460 only issue (checked with Kalyn, some of these realizations just don’t have future years for these experiments). Improved things a lot but I still have to look into the ones that are left.

# 
# for (draw in 1:Ndraws){
#   # The matching to the archive without ssp245 results
#   match_neighborhood(target_data = target_data,
#                      archive_data = archive_subset1,
#                      tol=0.1) %>%
#     # Convert them to a sample of individual Tgav Recipes:
#     permute_stitching_recipes(N_matches = 2000, matched_data = .,
#                               archive = archive_subset1)  %>%
#     # And stitch the recipes into tgavs:
#     stitch_global_mean(recipes=., data = tgav_data) %>%
#     # eliminate any Tgav time series that had NA's stitched in
#     group_by(stitching_id) %>%
#     mutate(drop = any(is.na(value))) %>%  
#     ungroup %>%
#     filter(drop == FALSE) %>%
#     select(-drop) %>% 
#     # add some identifying info
#     mutate(tol = 0.1,
#            archive = 'withoutSSP245') ->
#     rslts_1_0p1
#   
#   # Check for collpase in these results
#     rslts_1_0p1 %>%
#     select(year, value, stitching_id) %>%
#     spread(year, value) ->
#     wide
#   
#     # Collapse would have occured if any of these trajectories go through the same value
#     # in the same year -> if the length of distinct values in any column is less than
#     # the number of stitching_id's
#     counts <- apply(wide, 2, length)
#     
#     print(paste('Archive without ssp245; tol = 0.1: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))
#   
#     
#   match_neighborhood(target_data = target_data,
#                      archive_data = archive_subset1,
#                      tol=0.25) %>%
#     # Convert them to a sample of individual Tgav Recipes:
#     permute_stitching_recipes(N_matches = 2000, matched_data = .,
#                               archive = archive_subset1)  %>%
#     # And stitch the recipes into tgavs:
#     stitch_global_mean(recipes=., data = tgav_data) %>%
#      # eliminate any Tgav time series that had NA's stitched in
#     group_by(stitching_id) %>%
#     mutate(drop = any(is.na(value))) %>%  
#     ungroup %>%
#     filter(drop == FALSE) %>%
#     select(-drop) %>% 
#     # add some identifying info
#     mutate(tol = 0.25,
#            archive = 'withoutSSP245') ->
#     rslts_1_0p25
#   
#     # Check for collpase in these results
#     rslts_1_0p25 %>%
#     select(year, value, stitching_id) %>%
#     spread(year, value) ->
#     wide
#   
#     # Collapse would have occured if any of these trajectories go through the same value
#     # in the same year -> if the length of distinct values in any column is less than
#     # the number of stitching_id's
#     counts <- apply(wide, 2, length)
#     
#     print(paste('Archive without ssp245; tol = 0.25: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))
#     
#   
#   match_neighborhood(target_data = target_data,
#                      archive_data = archive_subset1,
#                      tol=0.5) %>%
#     # Convert them to a sample of individual Tgav Recipes:
#     permute_stitching_recipes(N_matches = 2000, matched_data = .,
#                               archive = archive_subset1)  %>%
#     # And stitch the recipes into tgavs:
#     stitch_global_mean(recipes=., data = tgav_data) %>%
#      # eliminate any Tgav time series that had NA's stitched in
#     group_by(stitching_id) %>%
#     mutate(drop = any(is.na(value))) %>%  
#     ungroup %>%
#     filter(drop == FALSE) %>%
#     select(-drop) %>% 
#     # add some identifying info
#     mutate(tol = 0.5,
#            archive = 'withoutSSP245') ->
#     rslts_1_0p5
#   
#   # Check for collpase in these results
#     rslts_1_0p5 %>%
#     select(year, value, stitching_id) %>%
#     spread(year, value) ->
#     wide
#   
#     # Collapse would have occured if any of these trajectories go through the same value
#     # in the same year -> if the length of distinct values in any column is less than
#     # the number of stitching_id's
#     counts <- apply(wide, 2, length)
#     
#     print(paste('Archive without ssp245; tol = 0.5: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))
#     
#   
#   # The matching to the archive with ssp245 results
#   match_neighborhood(target_data = target_data,
#                      archive_data = archive_subset2,
#                      tol=0.1) %>%
#     # Convert them to a sample of individual Tgav Recipes:
#     permute_stitching_recipes(N_matches = 2000, matched_data = .,
#                               archive = archive_subset2)  %>%
#     # And stitch the recipes into tgavs:
#     stitch_global_mean(recipes=., data = tgav_data) %>%
#      # eliminate any Tgav time series that had NA's stitched in
#     group_by(stitching_id) %>%
#     mutate(drop = any(is.na(value))) %>%  
#     ungroup %>%
#     filter(drop == FALSE) %>%
#     select(-drop) %>% 
#     # add some identifying info
#     mutate(tol = 0.1,
#            archive = 'withSSP245') ->
#     rslts_2_0p1
#   
#     # Check for collpase in these results
#     rslts_2_0p1 %>%
#     select(year, value, stitching_id) %>%
#     spread(year, value) ->
#     wide
#   
#     # Collapse would have occured if any of these trajectories go through the same value
#     # in the same year -> if the length of distinct values in any column is less than
#     # the number of stitching_id's
#     counts <- apply(wide, 2, length)
#     
#     print(paste('Archive with ssp245; tol = 0.1: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))
#     
#   
#   match_neighborhood(target_data = target_data,
#                      archive_data = archive_subset2,
#                      tol=0.25) %>%
#     # Convert them to a sample of individual Tgav Recipes:
#     permute_stitching_recipes(N_matches = 200, matched_data = .,
#                               archive = archive_subset2)  %>%
#     # And stitch the recipes into tgavs:
#     stitch_global_mean(recipes=., data = tgav_data) %>%
#      # eliminate any Tgav time series that had NA's stitched in
#     group_by(stitching_id) %>%
#     mutate(drop = any(is.na(value))) %>%  
#     ungroup %>%
#     filter(drop == FALSE) %>%
#     select(-drop) %>% 
#     # add some identifying info
#     mutate(tol = 0.25,
#            archive = 'withSSP245') ->
#     rslts_2_0p25
#   
#   # Check for collpase in these results
#     rslts_2_0p25 %>%
#     select(year, value, stitching_id) %>%
#     spread(year, value) ->
#     wide
#   
#     # Collapse would have occured if any of these trajectories go through the same value
#     # in the same year -> if the length of distinct values in any column is less than
#     # the number of stitching_id's
#     counts <- apply(wide, 2, length)
#     
#     print(paste('Archive with ssp245; tol = 0.25: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))
#   
#   
#   match_neighborhood(target_data = target_data,
#                      archive_data = archive_subset2,
#                      tol=0.5) %>%
#     # Convert them to a sample of individual Tgav Recipes:
#     permute_stitching_recipes(N_matches = 200, matched_data = .,
#                               archive = archive_subset2)  %>%
#     # And stitch the recipes into tgavs:
#     stitch_global_mean(recipes=., data = tgav_data) %>%
#      # eliminate any Tgav time series that had NA's stitched in
#     group_by(stitching_id) %>%
#     mutate(drop = any(is.na(value))) %>%  
#     ungroup %>%
#     filter(drop == FALSE) %>%
#     select(-drop) %>% 
#     # add some identifying info
#     mutate(tol = 0.5,
#            archive = 'withSSP245') ->
#     rslts_2_0p5
#   
#   # Check for collpase in these results
#     rslts_2_0p5 %>%
#     select(year, value, stitching_id) %>%
#     spread(year, value) ->
#     wide
#   
#     # Collapse would have occured if any of these trajectories go through the same value
#     # in the same year -> if the length of distinct values in any column is less than
#     # the number of stitching_id's
#     counts <- apply(wide, 2, length)
#     
#     print(paste('Archive with ssp245; tol = 0.5: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))
#   
#   
#   
#   write.csv(bind_rows(rslts_1_0p1, rslts_1_0p25, rslts_1_0p5,
#                       rslts_2_0p1, rslts_2_0p25, rslts_2_0p5) %>%
#               mutate(draw = draw),
#             paste0('generated_ensembles/generated_ensembles_ssp245_draw', draw, '.csv'), 
#             row.names = FALSE)
# }

sidenote - took about 10 min for each individual draw of these different generated ensemble.

Within each draw of a generated ensemble, we don’t necessarily have a ton of generated ensemble members.

Part of this is because of how we order and prioritize returning stitching recipes from a set of matches. For example, the most restrictive archive for matching, tol=0.1 and no ssp245 results. A couple ensemble members have a window with only one potential matched point in the neighborhood; Those get done first and if say realiation 12 and realization 13 both have the same archive point match for the same window, only one of them gets that point and the other doesn’t get any generated ensemble members. This is because we are enforcing the strongest possible restrictions to avoid envelope collapse. If we relaxed those, say by applying the pipeline to each target ensemble member separately, we’d get more matches but there would probably be some collapse acros target ensemble members:

match_neighborhood(target_data = target_data,
                   archive_data = archive_subset1,
                   tol=0.1) ->
  matched

kable(get_num_perms(matched)[[1]] )

Versus a still strict matching neighborhood (tol=0.1) but being able to draw from the ssp245 runs for the matching as well:

match_neighborhood(target_data = target_data,
                   archive_data = archive_subset2,
                   tol=0.1) ->
  matched2

kable(get_num_perms(matched2)[[1]] )

4 Validate - ensembles

Comparing the generated ensembles to the training ensemble with error metrics along the lines of Tebaldi et al 2020.

That is, we need the means and variances of both the target ensemble and the generated ensemble. Specifically, mean across time and ensemble members, variance across time and ensemble members.

# values of the target data
mean_target <- mean(original_data$value)
mean_target_hist <- mean(filter(original_data, year <= 2014)$value)
mean_target_future <- mean(filter(original_data, year > 2014)$value)
  
var_target <- sd(original_data$value)
var_target_hist <- sd(filter(original_data, year <= 2014)$value)
var_target_future <- sd(filter(original_data, year > 2014)$value)
  

# generated ensemble files
# read in all generated ensembles.
files <- list.files(path = 'generated_ensembles', pattern = 'generated_ensembles', full.names = TRUE) 

for (i in 1:length(files)){

    # Metrics for this draw #################################################################
  generated_ensemble_draw <- read.csv(files[i], stringsAsFactors = FALSE)
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'hist+future',
           mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  generated_ensemble_draw %>% 
    filter(year <= 2014) %>%
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'hist',
           mean_target = mean_target_hist,
           var_target = var_target_hist,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics_hist
  
  
  generated_ensemble_draw%>% 
    filter(year > 2014) %>%
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(value), 
              var_gen = sd(value)) %>%
    ungroup %>%
    mutate(time_period = 'future',
           mean_target = mean_target_future,
           var_target = var_target_future,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics_future
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup ->
    tmp
  
  bind_rows(tmp %>%
              left_join(metrics, by = c('tol', 'archive', 'draw')),
            tmp %>%
              left_join(metrics_hist, by = c('tol', 'archive', 'draw')) ,
            tmp %>%
              left_join(metrics_future, by = c('tol', 'archive', 'draw')))  ->
    plotting_tbl
  rm(tmp)
  
   write.csv(plotting_tbl %>%
              mutate(draw = draw),
            paste0('generated_ensembles/metrics_ssp245_draw', i, '.csv'), 
            row.names = FALSE)
}

# read in metrics for every draw
files <- list.files(path = 'generated_ensembles', pattern = 'metrics', full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive, time_period) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)

ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  facet_wrap(~time_period, nrow =3) +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0)

5 Validate - window jumps

Pull off the years that end each window and that start each window; these are going to be years we pay particular attention to.

window_start_years <- unique(archive_subset1$start_yr)
window_end_years <- unique(archive_subset1$end_yr)

We are particularly concerned about the ‘jumps’ at each of these years being ‘wrong’. And the idea is that they would be wrong relative to both the target data and, I think, the generated ensembles in between the jump years.

So let’s work with the jumps themselves for validating. Take the year to year diff within each realization (target or generated ensemble member). We’ll look at the diff and abs(diff). We’ll visualize the jumps with a scatter view and then get into more quantitative metrics. I spent a while going back and forth between starting from fundamentals like this and tracking down a package that looks at more sophisticated methods. As I started reading through packages, I ended up deciding to stick to fundamentals because they’ll translate from R to python better and will provide a more solid understanding of what’s going on before we look at anything more sophisticated.

# read in all generated ensembles.
files <- list.files(path = 'generated_ensembles', pattern = 'generated_ensembles', full.names = TRUE) 

generated_tgavs <- list()
for (i in 1:length(files)){
  generated_tgavs[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_tgavs <- do.call(rbind, generated_tgavs)


# make a helper function for calculating all the jumps because you have Nyears-1 jumps and
# dply functions don't like that. 
get_jumps <- function(data){
  data.frame(left_year = data$year[1:(length(data$year)-1)],
             jump = diff(data$value),
             abs_jump = abs(diff(data$value))) %>%
    mutate(variable = unique(data$variable),
           stitching_id = unique(data$stitching_id),
            tol = unique(data$tol),
           archive = unique(data$archive),
           draw = unique(data$draw),
           jump_year = if_else(left_year %in% window_end_years, 'window_transition', 'in_window')) ->
    x
return(x)
}

# calculate the yearly jumps
generated_list <- split(generated_tgavs, 
                        list(generated_tgavs$tol,
                        generated_tgavs$archive,
                        generated_tgavs$draw, 
                        generated_tgavs$stitching_id),
                        drop=TRUE)

lapply(generated_list, get_jumps) %>%
  do.call(rbind, .) ->
  generated_jumps

target_list <- split(original_data,
                     list(original_data$model, 
                          original_data$experiment,
                          original_data$ensemble),
                     drop = TRUE)

lapply(target_list, get_jumps) %>% 
  do.call(rbind, .) ->
  target_jumps

ggplot() + 
  geom_point(target_jumps, mapping = aes(x = left_year, y = jump), alpha = 0.5) +
  geom_point(filter(generated_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year, y = jump), alpha = 0.1, color='red') +
    facet_grid(tol~archive)

ggplot() + 
  geom_point(filter(target_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year, y = jump), alpha = 0.5) +
  geom_point(filter(generated_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year-1, y = jump), alpha = 0.1, color='red' ) +
    facet_grid(archive~tol) +
  ggtitle('note the years are staggered just for visual clarity')

# focusing on those transition years, what's the max value of the 
# target data for each of those years
filter(target_jumps,
       jump_year == 'window_transition') %>%
  group_by(left_year, variable, jump_year) %>%
  summarize(maxJump = max(abs_jump)) %>%
  ungroup ->
  target_max_values


ggplot() + 
  geom_point(filter(target_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year, y = abs_jump), alpha = 0.5) +
  geom_point(filter(generated_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year-1, y = abs_jump), alpha = 0.1, color='red' ) +
  geom_point(target_max_values, mapping = aes(x=left_year, y = maxJump), color = 'blue', shape = 4, size = 2) +
  geom_point(target_max_values, mapping = aes(x=left_year, y = maxJump + 0.44), color = 'blue', shape = 4, size = 2) +
  geom_point(target_max_values, mapping = aes(x=left_year, y = maxJump + 0.81), color = 'blue', shape = 4, size = 2) +
  geom_point(target_max_values, mapping = aes(x=left_year, y = maxJump + 1.54), color = 'blue', shape = 4, size = 2) +
    facet_grid(archive~tol) +
  ggtitle('note the years are staggered just for visual clarity. \n blue X = max target value in that year plus historic, future, total SDs as increase')

ggplot() + 
  geom_point(filter(target_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year, y = abs_jump), alpha = 0.5) +
  geom_point(filter(generated_jumps,
                    jump_year == 'window_transition'), mapping = aes(x = left_year-1, y = abs_jump), alpha = 0.1, color='red' ) +
  geom_hline(yintercept = max(target_jumps$abs_jump)) +
  geom_hline(yintercept = max(target_jumps$abs_jump + 0.44)) +
  geom_hline(yintercept = max(target_jumps$abs_jump + 0.81)) +
  geom_hline(yintercept = max(target_jumps$abs_jump + 1.54)) +
  facet_grid(archive~tol) +
  ggtitle('note the years are staggered just for visual clarity')

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))

focus_draw <- 1

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

focus_draw <- 2

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

focus_draw <- 27

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

ggplot() +
  geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'),  mapping = aes(x = jump, fill = id), alpha =0.2) +
  geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
  facet_grid(archive~tol) +
  ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))

# values of the target data
mean_target <- mean(target_jumps$jump)

var_target <- sd(target_jumps$jump)

  
for (drawID in 1:Ndraws){
  # Metrics for this draw #################################################################
  generated_jumps %>%
    filter(draw == drawID) ->
    generated_ensemble_draw 
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(jump), 
              var_gen = sd(jump)) %>%
    ungroup %>%
    mutate(mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup  %>%
    left_join(metrics, by = c('tol', 'archive', 'draw')) ->
    plotting_tbl
  
   write.csv(plotting_tbl %>%
              mutate(draw = drawID),
            paste0('generated_ensembles/jump_errors_allyears_ssp245_draw', drawID, '.csv'), 
            row.names = FALSE)
}


# Plot the first draw's errors:
plotting_tbl <- read.csv(paste0('generated_ensembles/jump_errors_allyears_ssp245_draw1.csv'),  
                         stringsAsFactors = FALSE)

kable(plotting_tbl)

ggplot(data = plotting_tbl, aes(x = E1, y = E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle('Error metrics for single Draw 1, comparing all \nyears of generated jumps to all years of target jumps')

# read in metrics for every draw
files <- list.files(path = 'generated_ensembles', pattern = 'jump_errors_allyears', full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)

ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle(paste('Aggregate errors (jumps for all years), aggregated across Ndraws ', unique(aggregate_metrics$total_draws) ))

# values of the target data
mean_target <- mean(target_jumps$jump)

var_target <- sd(target_jumps$jump)

  
for (drawID in 1:Ndraws){
  # Metrics for this draw #################################################################
  generated_jumps %>%
    filter(draw == drawID,
           jump_year=='window_transition') ->
    generated_ensemble_draw 
  
  # Values of the generated data and then the metrics
  generated_ensemble_draw %>% 
    group_by(tol, archive, draw) %>%
    summarize(mean_gen = mean(jump), 
              var_gen = sd(jump)) %>%
    ungroup %>%
    mutate(mean_target = mean_target,
           var_target = var_target,
           bias = abs(mean_target - mean_gen),
           E1 = bias/var_target,
           E2 = var_gen/var_target) ->
    metrics
  
  
  
  generated_ensemble_draw %>%
    select(stitching_id, tol, archive, draw) %>%
    distinct %>%
    group_by(tol, archive, draw) %>%
    summarise(n_stitched = n()) %>%
    ungroup  %>%
    left_join(metrics, by = c('tol', 'archive', 'draw')) ->
    plotting_tbl
  
   write.csv(plotting_tbl %>%
              mutate(draw = drawID),
            paste0('generated_ensembles/jump_errors_transitionyears_ssp245_draw', drawID, '.csv'), 
            row.names = FALSE)
}


# Plot the first draw's errors:
plotting_tbl <- read.csv(paste0('generated_ensembles/jump_errors_transitionyears_ssp245_draw1.csv'),  
                         stringsAsFactors = FALSE)

kable(plotting_tbl)

ggplot(data = plotting_tbl, aes(x = E1, y = E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle('Error metrics for single Draw 1, comparing window transition years of \ngenerated jumps to all years of target jumps')

# read in metrics for every draw
files <- list.files(path = 'generated_ensembles', pattern = 'jump_errors_transitionyears', full.names = TRUE) 

generated_metrics <- list()
for (i in 1:length(files)){
  generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)

generated_metrics %>% 
  group_by(tol, archive) %>%
  summarize(aggregate_E1 = mean(E1),
            aggregate_E2 = mean(E2)) %>%
  ungroup %>%
  mutate(total_draws = length(unique(generated_metrics$draw))) ->
  aggregate_metrics

kable(aggregate_metrics)

ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
  geom_point() +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 0) + 
  ggtitle(paste('Aggregate errors comparing window transition years of \ngenerated jumps to all years of target jumps \naggregated across Ndraws ', unique(aggregate_metrics$total_draws) ))

6 Next?

Need to look at other scenarios and eventually more models. The Tgav analysis for sure seems to be the way to give us hyper parameters but I can’t tell yet where the rule of thumb will come from. Or if the best hyperparamters will be the same across scenarios and/or models - want the largest possible tol that still gives good emulation in all uses, basically.

• look at the NA thing; either there’s more collapse free ensemble members we could be getting or there’s a test I can/should add into the function. Seems to be largely coming in around the 2014-2015 ssp434 ssp460 thing. -> KD and I checked briefly and some of the realizations for these experiments don’t have future data at all ; That’s been fixed but there’s a few NAs getting introduced in matching left. Not fatal, and not surprising I missed it since I didn’t use those two scenarios a ton during my development. I’ll be tracking it down.

• parallelize and test locally, then launch on pic to get actual 1000s of draws for validating the method. Should be pretty easy with a dopar over the draw id value, everyone should be working separately and saving things with distinct ids that go with them when the job is dispatched to a node. Can include the metrics on tgav and year to year diffs in the dispatch so that all we have to do is read everything in to aggregate and plot. Also functionalize the errors to speed up

• Look at SSP585 specifically

• maybe something more sophisticated than just comparing distributions of the jumps. Maybe look at the relative frequency of how often a generated jump is outside the range of target jumps. We’d need a large number of draws for that,and it isn’t even particularly sophisticated.

7 What about monthly, gridded data?

Assuming we are happy with the above style of analysis looking at both the general method behavior (many draws of full generated ensembles) + individual use case behavior (one draw of full generated ensemble).

• Doing the method test on gridded data might be too challenging. 1000s of annual Tgav time series is very different from all the opening and closing of the netcdfs to construct monthly (or daily) gridded time series. We would at least not have to repeat the hyperparameter investigation on the gridded data though, which is one of the reasons the method behavior validation was desirable. I think it’s reasonable to say ‘we’ve set the hyperparamters based on validating Tgavs’. Could get maybe 100 draws on the gridded data under just the final set of hyperparameters then. Definitely not enough out of the (total possible permutations choose Nmin collapse free trajectories) solution space but we can cross that bridge when we have a better sense of emulation times after a little bit of optimizing maybe? Have to consider that we’ll have to repeat analysis for all RCPS, multiple ESMs.

• from the gridded perspective, if it was annual data, I think we could look at the metrics on the actual values and on the year to year jump values in each grid cell and have summaries of those grid cell behaviors to know if we’re happy or not, basically repeating what we did above and either doing an area weighted average to get global E1, E2 values or just look at the quantiles of the E1, E2 values over grid cells. Should be good enough to cover time and we could look at EOFs and power spectra to maybe cover space (have to think a lot more about that).

• KS test might be good to do by grid cell on the distributions of jumps I suppose, to supplement the temporal analysis.

• I don’t know how to think about the monthly factor though. Would just looking at the month to month jumps be enough on the temporal side? And focus on maybe the Dec to Jan jump where we paste? Or the Jan-Dec and then Jan-Dec at a jump? Spatially, I’m not sure. I guess EOFs and power spectra could still be done, but I need to read more about what they really mean for monthly data instead of annual.